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We present the first complete calculation of the parameter / NL , a quantity introduced to char- 
acterize the extent of non-Gaussianity, for a variety of single field inflationary models that lead to 
features in the scalar power spectrum. The calculation is performed numerically by means of a new, 
efficient and accurate Fortran code that can evaluate all the contributions to the bi-spectrum in any 
configuration. We consider different sets of models that lead to similar features in the scalar power 
spectrum, and investigate if (viz. / NL evaluated in the equilateral configuration) can help us 
discriminate between the models. We find that certain differences in the background dynamics — 
reflected in the behavior of the slow roll parameters — can lead to a reasonably large difference in 
the generated by the models. We close with a discussion on the implications of the results we 
obtain. 

PACS numbers: 98.80.Cq, 98.70. Vc 
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Primordial spectra with features: Most single field infla- 
tionary models naturally lead to an extended period of 
slow roll, and hence to nearly scale invariant primordial 
spectra, which seem to be fairly consistent with the re- 
cent data on the anisotropics in the Cosmic Microwave 
Background (CMB) [1] as well as other observational con- 
straints. However, it has been repeatedly noticed that 
certain features in the inflationary scalar power spec- 
trum can improve the fit to the CMB data at the cost of 
some additional parameters (see, for instance, the recent 
work [2] and the long list of references cited therein). 
Though the statistical significance of these features re- 
main to be understood satisfactorily (see, for example, 
Refs. [3]), they gain importance from the phenomenolog- 
ical perspective of comparing the models with the data, 
because only by a smaller class of single field inflationary 
models, which allow for departures from slow roll, can 
generate them. 

Over the last half-a-dozen years, it has been increas- 
ingly realized that the detection of non-Gaussianities in 
the primordial perturbations can considerably help in 
constraining the inflationary models (see Refs. [4, 5]; for 
early efforts in this direction, see Refs. [6]). In particu- 
lar, the detection of a high value for the / NL parameter 
that is used to describe the extent of non-Gaussianity 
can rule out a wide class of models. If, indeed, the ex- 
tent of non-Gaussianity proves to be as large as the mean 
values of / NL arrived at from the recent WMAP data [7- 
9], then canonical scalar field models that lead to slow 
roll inflation and nearly scale invariant primordial spec- 
tra will cease to be consistent with the data. But, inter- 
estingly, demanding the presence of features in the scalar 
power spectrum seems to gcnerically lead to larger non- 
Gaussianities (see, for example, Refs. [10]). Therefore, 
features may offer the only route (unless one works with 
non- vacuum initial states [11]) for the canonical scalar 



fields to remain viable if / NL turns out to be significant. 

The above discussion raises two important issues. 
Firstly, if indeed the presence of features turns out 
to be the correct reason behind possibly large non- 
Gaussianities, can we observationally identify the correct 
underlying inflationary scenario, in particular, given the 
fact that different models can lead to similar features 
in the scalar power spectrum? In other words, to what 
extent can the non-Gaussianity parameter / NL help us 
discriminate between the inflationary models that per- 
mit features? To address this question, we shall con- 
sider a few typical inflationary models leading to fea- 
tures, assuming that they can be viewed as representa- 
tives of such a class of scenarios. Concretely, we shall 
consider the Starobinsky model [12] and the punctuated 
inflationary scenario [13], both of which result in a sharp 
drop in power at large scales that is followed by oscilla- 
tions. We shall also study large and small field models 
with an additional step introduced in the inflaton poten- 
tial [10, 14]. The step leads to a burst of oscillations in 
the scalar power spectrum which improve the fit to the 
outliers near the multipole moments of I = 22 and 40 in 
the CMB anisotropies. We shall also consider oscillating 
inflaton potentials such as the one that arises in the ax- 
ion monodromy model which lead to modulations in the 
power spectrum over a wide range of scales [2, 10, 15]. 
In Fig. 1, we have illustrated the scalar power spectrum 
that arises in these different models. 

The second issue pertains to the calculation of non- 
Gaussianities in models where the slow roll approxima- 
tion is not satisfied. Usually, the slow roll approximation 
is utilized to arrive at analytical expressions for the non- 
Gaussianity parameter / NL . Clearly, this is no longer pos- 
sible when departures from slow roll occur. In this work, 
we shall use a new Fortran numerical code to evaluate 
the non-Gaussianities in such situations. Although, some 
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FIG. 1: The scalar power spectrum in the various inflationary 
models that we consider. The parameters of the Starobin- 
sky model [12] has been chosen such that the resulting power 
spectrum closely resembles the spectrum that arises in the 
punctuated inflationary scenario which is known to lead to 
an improved fit to the CMB data [13]. While the models with 
a step [14] lead to a burst of oscillations over a specific range 
of scales, inflaton potentials with oscillating terms produce 
modulations over a wide range of scales in the power spec- 
trum [2]. The inset highlights the differences in the various 
power spectra over a smaller range of scales. 

partial numerical results have already been published in 
the literature, we believe that it is for the first time that a 
general (we shall restrict ourselves to the equilateral case 
in this letter, but the code can compute for any configu- 
ration), and efficient (that can arrive at results within a 
few minutes) code has been put together. Moreover, as 
we shall demonstrate, the code can also compute all the 
different contributions to the bi-spectrum. 

The plan of this letter is as follows. We shall first 
quickly describe the essential details pertaining to the 
evaluation of the non-Gaussianity parameter / NL in in- 
flationary models involving a single canonical scalar field. 
We shall then describe the method that we adopt to nu- 
merically compute the parameter / NL in the equilateral 
limit, and illustrate the extent of accuracy of the com- 
putations by comparing our numerical results with the 
analytical results that have recently been obtained in the 
case of the Starobinsky model [16]. Subsequently, we 
shall present the main results, and compare the / NL that 
arise in the various models of our interest. We finally 
conclude with a brief discussion on the implications of 
our results. 

/ NL — Essentials: Let us quickly recall the essentials. 
The scalar power spectrum V s (k) and the bi-spectrum 
6 s (ki, k 2 , k 3 ) are denned in terms of the two and three 
point correlation functions of the Fourier modes of the 
curvature perturbation 1Z as follows [7]: 

(£k£ P > - ^P B (*)« (3) (k + p), (1) 
(K^K^K^) = (2tt) 3 ^(k^k^ka) 

X5 (3) (k!+k 2 + k3), (2) 

where k = |k|. The dimcnsionless non-Gaussianity pa- 



rameter / NL is introduced through the relation 

U = K G - (K G ) 2 , (3) 
5 

where 7Z G denotes the Gaussian contribution to 7Z. Upon 
using this relation, and the definitions of the power spec- 
trum and the bi-spectrum above, one can arrive at the 
expression for the parameter / NL in terms of the bi- 
spectrum and the power spectrum. In the equilateral 
limit of our interest, i.e. when hi — k 2 = k 3 = k, it is 
found to be 

_1Q _J_ fc 6 G(fc) 
/nl 9 (2tt) 4 Vl(k) ' 1 ' 

where G(k) = (2 7r) 9 / 2 B s (k). It ought to be stressed here 
that the non-Gaussianity parameter / NL is usually intro- 
duced through Eq. (3) with the local limit in mind. The 
above expression for f^3 has been analogously defined in 
the equilateral limit. 

There now exists a standard procedure for evaluating 
the scalar bi-spectrum in inflationary models [4, 5, 10]. 
The quantity G(k), evaluated towards the end of inflation 
at the conformal time, say, r\ = rj e , can be written as [16] 

6 

G(k) ^ m* [fk faO G c (k) + a: 3 (77c) g*(k)] 

C=l 

+ G 7 (k), (5) 

where the quantities Q c (k) are integrals that correspond 
to six terms that arise in the action at the third order 
in the perturbations [4, 10], while fk are the modes as- 
sociated with the curvature perturbation. For instance, 
the fourth term Qi{k), which proves to be the dominant 
contribution to / NL in certain situations when deviations 
from slow roll occur, is described by the following integral 
over conformal time: 

e 4 (fc)=37 f° dr^e^/f rC (6) 

where a denotes the scale factor, t\ and e 2 are the first 
two slow roll parameters, the overprime represents differ- 
entiation with respect to the conformal time coordinate, 
while rji is an initial time when, say, the largest mode of 
interest is well inside the Hubble radius. The additional, 
seventh term arises due to a field redefinition [4, 10]. 
The numerical computation of /^j : We shall now briefly 
outline the methods that we adopt to numerically evolve 
the equations governing the background and the pertur- 
bations, and eventually evaluate the inflationary scalar 
power and bi-spectra. 

We solve the background as well as the perturbation 
equations using a Bulirsch-Stoer algorithm with an adap- 
tive step size control routine [17]. We shall treat the 
number of e-folds as the independent variable, which al- 
lows for a more efficient and accurate computation. Since 
we shall be focusing on the equilateral limit of the bi- 
spectrum, we can evolve each of the modes of interest 
independently and compute the resulting f°^. We im- 
pose the standard Bunch-Davies initial conditions on the 
perturbations when the modes are well inside the Hub- 
ble radius H [say, when (k/aH) = 100], and evolve them 
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FIG. 2: The quantities k times the absolute values of G1+G3 
(in green), G2 (in red), G4 (in blue), G5 + Ge (in purple) and 
G7 (in dark green) have been plotted as a function of k/ko 
for the Starobinsky model. Note that ko is the wavenumber 
which leaves the Hubble radius when the scalar field crosses 
the break in the potential. The solid curves represent the ana- 
lytical expressions that have been obtained recently [16], while 
the dashed curves denote the numerical results computed us- 
ing our Fortran code. The dots of an alternate color denote 
the corresponding numerical values that have been arrived at 
independently using a Mathematica [19] code. We find that 
the numerical results match the analytical results exception- 
ally well in the case of the crucial, dominant contribution to 
the / NL , viz. due to G4. 

until very late times 1 . The power spectra displayed in 
Fig. 1 have been evaluated at super Hubble scales, when 
the amplitude of the curvature perturbations have frozen 
in [typically, when (k/aH) ~ 10~ 5 ]. We carry out the 
integrals involved in arriving at the bi-spectrum using 
the method of adaptive quadrature [18]. We compute 
the integrals from the time when the initial conditions 
are imposed on the modes till the time when they are 
well outside the Hubble radius. The integrals Q n contain 
a cut off in the sub Hubble limit, which singles out the 
perturbative vacuum. Numerically, the presence of such 
a cut off is fortunate since it controls the contributions 
due to the continuing oscillations that would otherwise 
occur. Generalizing the cut off that is often introduced 
analytically in the slow roll case, we impose a cut off of 
the form exp — [S (k/aH)], where S is a small parame- 
ter. The various tests that we have carried out indicate 
that the results are fairly robust to changes in the value 
of S and the limits of the integration, provided we inte- 
grate from sufficiently deep inside the Hubble radius till 



This is so apart from the case of the axion monodromy model 
wherein the modes have to be integrated from a suitably early 
initial time so that the resonance that occurs in these models 
due to the oscillations in the potential is captured [10, 15]. 



suitably late times. 

Comparison with the analytical results: The Starobinsky 
model consists of a linear inflaton potential with an 
abrupt change in the slope [12]. The change in the slope 
causes a brief period of fast roll which leads to sharp 
features in the scalar power spectrum (as illustrated in 
Fig. 1). It was known that, for certain range of parame- 
ters, one could evaluate the scalar power spectrum ana- 
lytically in the Starobinsky model, which matches the ac- 
tual, numerically computed spectrum exceptionally well. 
Interestingly, two of us have recently shown that, in the 
equilateral limit, the model allows the analytic evalua- 
tion of the scalar bi-spectrum too [16]. Before we go 
on to consider other models, we shall compare the nu- 
merical results we obtain with the analytical results that 
are available in the case of the Starobinsky model. In the 
equilateral limit, the contributions due to G\ and G3 and 
G5 and Gq turn out to be of the same form. In Fig. 2, 
we have plotted the numerical as well as the analytical 
results for the functions G\ + G3, G2, G4, G5 + Gq and 
Gj for the Starobinsky model. We have plotted for pa- 
rameters of the model for which the analytical results 
are considered to be a good approximation [16]. It is 
evident from the figure that the numerical results match 
the analytical ones very well. Importantly, the agreement 
proves to be excellent in the case of the dominant contri- 
bution G4 and, as one would expect (since it involves no 
integrals), in the case of G7. 

Results: In Fig. 3, we have plotted the various contri- 
butions, viz. Gi + G3, G2, G4, G5 + Gq and G7 (in 
the equilateral limit) for punctuated inflation [13], the 
chaotic model with a step [14], and the axion monodromy 
model which contains oscillations in the inflaton poten- 
tial [2, 15]. These plots and the one in previous figure 
clearly point to the fact that, barring the case of punc- 
tuated inflation wherein G7 becomes dominant at large 
wavenumbers, it is the quantity G4 that contributes the 
most to the scalar bi-spectrum when deviations from slow 
roll occur [10]. In Fig. 4, we have plotted the due to 
the dominant contribution(s) that arises in the various 
models that we have considered. It is clear from this 
figure that, while in certain cases can prove to be a 
good discriminator, it cannot help in others, and its abil- 
ity to discriminate depends strongly on the differences in 
the background dynamics. For instance, the evolution of 
the first two slow roll parameters are very similar when 
a step is introduced in cither the quadratic potential or 
a small field model [14]. Hence, it is not surprising that 
the f oq behaves in a similar fashion in both these mod- 
els. Whereas, f®3 proves to be substantially different in 
punctuated inflation and the Starobinsky model. This 
difference can be attributed to the fact that, while in the 
Starobinsky model, the first slow roll parameter remains 
small throughout the evolution, it grows above unity for a 
very short period (leading to brief interruption of the ac- 
celerated expansion) in the punctuated inflationary see- 
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FIG. 3: The set of quantities k 6 \G n (k)\ plotted as in the pre- 
vious figure with the same choice of colors to represent the 
different G n (k). The figures on top, in the middle and at the 
bottom correspond to punctuated inflation, the quadratic po- 
tential with a step and the axion monodromy model, respec- 
tively, and they have been plotted for values of the parameters 
that lead to the best fit to the WMAP data [13-15]. In the 
middle figure, the dashed lines correspond to the quadratic 
potential when the step is not present. 

nario. This departure from inflation results in a sharp 
drop in power, which in turn leads to a sharp rise in f™*, 
possibly indicating that punctuated inflation is already 
ruled out. Similarly, we find that is rather large in the 
axion monodromy model in contrast to the case wherein 
the conventional quadratic potential is modulated by an 
oscillatory term (we have not displayed the results for 
the latter for want of space). The large value of that 
arises in the monodromy model can be attributed to the 
resonant behavior encountered in the model [10, 15]. 

Discussion: In this work, we have been interested in ex- 
amining the power of the non-Gaussianity parameter / NL 
to discriminate between various single field inflationary 
models involving the canonical scalar field. With this 
goal in mind, using a new numerical code which can effi- 
ciently calculate the bi-spectrum for any triangular con- 
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FIG. 4: A plot of corresponding to the various models 
that we have considered. The figure at the top contains the 
absolute value of , plotted on a logarithmic scale (for con- 
venience in illustrating the extremely large values that arise), 
in the Starobinsky model and the punctuated inflationary sce- 
nario. The figure in the middle contains /° q for the cases 
wherein a step has been introduced in a quadratic potential 
and a small field model. The figure at the bottom corre- 
sponds to that of the axion monodromy model. In fact, we 
have also evaluated the for the case of quadratic poten- 
tial modulated by sinusoidal oscillations, which too leads to 
continuing, periodic features in the scalar power spectrum. 
However, we find that the in such a case proves to be 
rather small (of the order 10~ 3 or so). As we had mentioned 
before, these sets of models lead to scalar power spectra with 
certain common characteristics. Needless to say, while is 
considerably different in the first and the last sets of models, 
it is almost the same in the case of models with a step. These 
similarities and differences can be attributed completely to 
the background dynamics. 

figuration, we have evaluated the quantity f*** in a slew of 
models that generate features in the scalar perturbation 
spectrum. We find that the amplitude of proves to 
be rather different when the dynamics of the background 
turns out reasonably different, which, in retrospect, need 
not be surprising at all. For instance, models such as the 
punctuated inflationary scenario which lead to very sharp 
features in the power spectrum also lead to substantially 
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FIG. 5: The various contributions to the bi-spectrum in the 
case of the Starobinsky model for a different set of parameters 
than those that have been plotted earlier in Fig. 2. Note that, 
the hierarchy of the various contributions is different with G2 
now being of the same order as G4 at large wavenumbers. 



large / NL . Such possibilities can aid us discriminate be- 
tween the models to some extent. We had focused on 
evaluating the quantity / NL in the equilateral limit. It 
will be worthwhile to compute the corresponding values 
in the other limits, such as the squeezed and the orthog- 
onal limits as well. We are currently investigating such 
issues [20]. 

We would like to conclude by highlighting two points. 
The first point concerns the hierarchy of the various con- 
tributions to the bi-spectrum. It has often been said that 
it is the term G4 that leads to the dominant contribution 
when deviations from slow roll occur. But, this need 
not be the case. In a recent work [16], two of us had 
illustrated as to how the contribution due to the second 
term can be as large as the fourth term under certain 
conditions in the case of the Starobinsky model. We 
have since been able to numerically confirm this result 
with our new code, a result which we have illustrated 
above in Fig. 5. Further, note that G7 dominates G4 at 
large wavenumbers in punctuated inflation. These results 
point to the fact that, when departures from slow roll 
occur, it becomes imperative that one evaluates all the 
contributions to the bi-spectrum, as we have done here. 
Secondly, the next logical step would be to compute the 
corresponding CMB bi-spectrum, an issue which we have 
not touched upon as it is beyond the scope of the cur- 
rent work. While tools seem to be available to evaluate 
the CMB bi-spectrum based on the first order brightness 
function, the contribution due to the brightness function 
at the second order remains to be understood satisfacto- 
rily (in this context, see Ref. [21] and the last reference 
in Ref. [8]). These seem important aspects worth inves- 
tigating closer. 
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